Pleiotropy promotes the evolution of inducible immune responses in a model of host-pathogen coevolution

Components of immune systems face significant selective pressure to efficiently use organismal resources, mitigate infection, and resist parasitic manipulation. A theoretically optimal immune defense balances investment in constitutive and inducible immune components depending on the kinds of parasites encountered, but genetic and dynamic constraints can force deviation away from theoretical optima. One such potential constraint is pleiotropy, the phenomenon where a single gene affects multiple phenotypes. Although pleiotropy can prevent or dramatically slow adaptive evolution, it is prevalent in the signaling networks that compose metazoan immune systems. We hypothesized that pleiotropy is maintained in immune signaling networks despite slowed adaptive evolution because it provides some other advantage, such as forcing network evolution to compensate in ways that increase host fitness during infection. To study the effects of pleiotropy on the evolution of immune signaling networks, we used an agent-based modeling approach to evolve a population of host immune systems infected by simultaneously co-evolving parasites. Four kinds of pleiotropic restrictions on evolvability were incorporated into the networks, and their evolutionary outcomes were compared to, and competed against, non-pleiotropic networks. As the networks evolved, we tracked several metrics of immune network complexity, relative investment in inducible and constitutive defenses, and features associated with the winners and losers of competitive simulations. Our results suggest non-pleiotropic networks evolve to deploy highly constitutive immune responses regardless of parasite prevalence, but some implementations of pleiotropy favor the evolution of highly inducible immunity. These inducible pleiotropic networks are no less fit than non-pleiotropic networks and can out-compete non-pleiotropic networks in competitive simulations. These provide a theoretical explanation for the prevalence of pleiotropic genes in immune systems and highlight a mechanism that could facilitate the evolution of inducible immune responses.

Introduction Organismal immune systems have been honed over evolutionary time to balance the benefits of infection resistance against the energetic and pathological costs of mounting and maintaining an immune response [1]. Optimal investment in immunity varies based on organismal life histories [2][3][4], but often features constitutive immune activity [5,6], wherein a defense mechanism is always active even when no threat is present. In humans, such quiescent immune activity can exceed 15% of basal metabolic demands, making the immune system one of the most energetically costly body systems, even in the absence of pathogens [7]. In contrast to constitutive responses, inducible immune responses are inactive until a threat is detected, at which point they quickly come online, clear the threat, and undergo rapid suppression by negative regulators [8]. As such, inducible immunity is a smaller drain on host fitness than constitutive immunity, but the time needed to activate the response may cede the temporal advantage against rapidly replicating pathogens [6,9].
Previous studies attempting to define optimal investment in constitutive and inducible immune responses have suggested that high parasite virulence and predictability, including frequent exposure and fixed within-host growth rates [10], favor the evolution of constitutive immunity [6]. Inducible immunity is favored when infection rates are low, parasite growth rates are varied [10], immune effector proteins can be produced rapidly [6], and when parasites are not co-evolving with the host [11]. If a parasite is coevolving with a host genotype, however, the host may have to contend with parasite-mediated manipulation, including suppressed inflammation, receptor interference, and degraded antigen presentation and processing [12]. A recent study on the evolution of immune signaling networks suggested that selection may favor constitutive defenses even when an inducible response is theoretically optimal because inducible architectures and their evolutionary intermediates may be less robust to this parasitic manipulation [11]. By identifying structural factors that promote the evolution of theoretically suboptimal network architectures, this study also provides a framework to investigate other genetic and genomic factors that could facilitate or hinder the evolution of inducible immune responses.
A ubiquitous but puzzling property of many immune signaling networks, for example, may either reinforce or resolve these constraints on the evolution of inducibility. Genetic pleiotropy, the phenomenon in which a single gene affects multiple distinct phenotypes, has been described in immune signaling networks from a variety of taxa. In plants, for example, DELLA proteins regulate network crosstalk between gibberellin-mediated growth and jasmonic-acid mediated defense, affecting two traits that compete for limited energetic resources [13]. In insects, the Toll pathway is crucial to both developmental patterning in fly embryos and defense against bacterial and fungal infections [14]. Recent work has underscored the high frequency of pleiotropy in insect immune and developmental signaling pathways [4,15] and a study of the human genome identified an overabundance of pleiotropy in the immune system [16]. Thus, pleiotropy appears to be a common property of immune signaling networks, but it is not clear whether and when it might constrain or benefit adaptive evolution.
Pleiotropy could exacerbate antagonism among traits at both within-host and evolutionary scales, and many empirical lines of evidence emphasize the potential for constraint. For example, the plasticity of pleiotropic gene expression in response to environmental perturbations is restricted relative to non-pleiotropic counterparts [17]. Pleiotropy can also impede the efficacy of selection, as purifying selection to retain one phenotype can interfere with substitutions that might provide an adaptive benefit for another [15,18,19]; this effect is magnified between traits that contribute equally to organismal fitness [20]. We further hypothesize that the elevated rate of purifying selection on some pleiotropic genes [15,21] and their abundance in the immune system could limit immune network adaptability and present stable avenues of manipulation to parasites.
Pleiotropy is not always deleterious, however, as computational [22] and genomics studies [23,24] of trait evolution demonstrate that pleiotropy commonly arises as a neutral process associated with the generation of new traits. Pleiotropic compression of genomes could also be vital for genome size management, which has implications for organismal fitness [19,25]. While these aspects of pleiotropy positively affect organismal fitness, it remains unclear what other beneficial changes to network evolution may accompany the incorporation of pleiotropic signaling proteins. To address this open question, we constructed an agent-based model of host parasite co-evolution featuring immune signaling networks and parasites that manipulate host signaling.
Our networks include three fundamental components: detectors, signaling proteins, and effectors. In an immunological context, detectors are pattern recognition receptor proteins like peptidoglycan recognition proteins and Toll-like receptors that directly sense parasites and initiate a response. Signaling proteins can be part of extracellular or intracellular signaling cascades that relay information from detector proteins to effector proteins. Effector proteins act to kill parasites or inhibit growth [26]; while some effectors across the tree of life exhibit behavior that is exclusively constitutive or inducible, here we focus on effectors like antimicrobial peptides (AMPs) in insects that can be both constitutively produced and induced through the transduction of a signal through a network after detection of a parasite [27]. We represent the proteins and protein-protein interactions that compose a signaling network as a collection of nodes and edges respectively. The pleiotropic protein explicitly contributes to immune signaling but also implicitly regulates another complex trait like development [28,29], such that its evolution is inherently constrained. Here we present four complementary implementations of this pleiotropic protein: 1) random connections between the pleiotropic protein and other proteins in the network that are fixed across evolutionary time; 2&3) a single connection from the pleiotropic protein that up-or downregulates the effector as a side effect of its function in development; and 4) the protein is allowed to evolve but under reduced mutation rates relative to other network components, capturing reduced evolutionary rates associated with purifying selection (Fig 1).
We first investigated whether these implementations of pleiotropy favor distinct patterns of constitutive and inducible immunity over evolutionary time relative to non-pleiotropic networks, keeping in mind that the predictability of infection could influence the outcome. To understand why pleiotropy might be such a common property of immune systems, we then investigated whether pleiotropic networks are capable of outcompeting non-pleiotropic signaling networks over evolutionary time and defined the network properties associated with successful networks. Our results provide a new framework for understanding the evolutionary maintenance of pleiotropy in immune systems and could be applied broadly to understand the evolution of immune signaling networks or the evolution of pleiotropic network architecture among the myriad traits that influence organismal fitness.

Results
Our agent-based model of signaling network co-evolution features hosts, defined by their immune networks, and parasites that disrupt host signaling to improve their own reproductive success. We present data from 100 simulations for each form of pleiotropic constraint and a non-pleiotropic case at three levels of lifetime infection risk [10%, 50%, 90%], with all simulations consisting of 500 generations. The initial host population was created with random connections between signaling proteins, and each host had the same implementation of pleiotropy within a scenario. Host fitness was evaluated based on the following components: pre-infection effector levels, capturing the cost of constitutive investment in immunity, post-infection effector levels, capturing the cost of excess investment in immunity, cumulative magnitude of parasite abundance, capturing the cost to host fitness for not managing their infection, and network size, capturing the fitness costs associated with a large genome (see Materials and Methods section 5). Parasite fitness was determined by the cumulative magnitude of infection over time within each host. Host evolution allowed for the duplication and deletion of signaling proteins as well as adding, removing, or altering connections between any two proteins in the network. Pleiotropic connections were immutable except in the case of slowed evolution, but paralogs created by the duplication of a pleiotropic protein were not similarly constrained. No direct connection between the detector and effector was allowed. Parasites were allowed to alter the identity and regulatory activity of their targeted signaling protein over evolutionary time. To evaluate the degree to which hosts could successfully mount an immune defense during infection (network robustness) we calculated the mean absolute difference in effector abundance between a host with an intact signaling network and that same host with a single signaling protein removed from the network. Both the intact and the knockout host where infected by a parasite that could grow in the host but could not manipulate host signaling. To compare the dynamics of the intact and knockout networks, we calculated the mean absolute difference between intact immune effector abundance and knockout immune effector abundance. To evaluate fundamental properties associated with these networks we measured the number of nodes in a network, total connectivity (the number of edges in the network divided by the total possible number of edges), and the distinct paths from the detector to effector protein, where a distinct path does not include a signaling protein used in any other path. We evaluated the fitness benefits of specific pleiotropic implementations versus the non-pleiotropic hosts by conducting competitive simulations where half of the initial population was pleiotropic, and half was non-pleiotropic. Hosts reproduced as normal, passing their pleiotropic status to their offspring, and a competitive simulation ended when one population was driven to extinction, or 1000 generations had passed ending the simulation in a draw.

A higher infection risk favors the evolution of stronger and inducible immune signaling
To determine how pleiotropic proteins may alter immune signaling network evolution, we first needed a baseline understanding of how non-pleiotropic networks evolve. We modeled the coevolution of non-pleiotropic signaling networks at three chances of infection to uncover the relationship between parasite prevalence and host immune evolution. When the individual chance of infection in a simulation was 10%, most hosts evolved predominantly constitutive immune responses (Fig 2A) with minimal investment in immunity ( Fig 2B). As the chance of infection increased, the likelihood of hosts developing a mixed-strategy immune response also increased, though predominantly inducible immune responses were still rare. Compared to the 10% chance of infection, hosts that evolved in the 50% chance of infection condition invested heavily in immunity, but this investment was still predominantly constitutive in nature. Hosts that evolved with a 90% chance of infection showed similar investment in immunity to those that evolved with a 50% chance of infection. When inducible responses evolved, they tended to have higher peak immune effector abundance than constitutive responses ( Fig 2B).
The 10% chance of infection produced hosts that were robust to parasitic interference, as defined by parasitic manipulation of host immune signaling to suppress immune effector abundance. There was a negligible difference, moreover, in immune effector dynamics between knockouts and intact networks during infection. As the chance of infection increased, we observed a larger discrepancy between intact and knockout networks (Fig 3). Network size and the distinct paths from the detector to the effector increased with the chance of infection, but network connectivity stayed consistent across all chances of infection (Fig K-M in S1 Text). Finally, inducibility was not favored when the chance of infection was low, with winners of competitions being no more inducible than losers (Fig 4).

Pleiotropy alters immune signaling network evolution
Broadly, the immune responses generated by pleiotropic networks were identical to those generated by non-pleiotropic networks, not just in the distribution of their immune response on the constitutive-inducible spectrum but also in the peak amount of immune effector deployed during infection (Fig 2A and 2B). When accounting for the differences in peak effector abundance between constitutive and inducible immune responses, there was no difference in the Normalized probability density function of the proportion of host immune responses that are induced by parasites. The x-axis shows the percent of the response that is induced by parasites, with the left-hand side being 0% of response induced, to 100% induced responses on the right. The y-axis corresponds to the relative likelihood of finding an immune response in the specified population that is X% induced. Non-pleiotropic networks are represented in green and pleiotropic networks in blue. All plots on the same row have the same chance of infection (10%, 50%, or 90% descending), and plots in the same column compare the non-pleiotropic network against the indicated pleiotropic constraint. Presented in each plot is the Pearson magnitude of immune effector abundance between pleiotropic and non-pleiotropic networks as a function of infection chance (Fig 2B). In four scenarios, however (fixed upregulation when the chance of infection was 10% or 90%, fixed downregulation when the chance of infection was 50% or 90%), the Pearson correlation coefficient between the non-pleiotropic immune response probability density function and the pleiotropic immune response probability density function showed little to no correlation (corr.�±.2), suggesting that these pleiotropic populations evolved different immune dynamics than those evolved by the non-pleiotropic populations.
When investigating active immune effector abundance changes associated with signaling protein knockouts, in most conditions the loss of the pleiotropic protein did not significantly alter the abundance of active immune effectors relative to the removal of any other signaling protein (Fig 3). Fixed downregulation is the only condition in which this was not true, indicating that theses hosts were reliant on the action of the pleiotropic signaling protein to produce their evolved immune response. On average, pleiotropic hosts were nearly as fit as non-pleiotropic hosts, with fixed downregulatory pleiotropy resulting in markedly higher absolute fitness than all other conditions when the chance of infection was 50% or 90% (

Pleiotropy can imbue competitive benefits to organisms
We used two classes of competitive simulations to evaluate the fitness of the pleiotropic populations relative to the non-pleiotropic ones. Unevolved competition began immediately following host initiation and evolved competition began after host populations had evolved for 250 generations in isolation (i.e., in the absence of a different type of pleiotropic implementation). In unevolved competition, we observed that the pleiotropic upregulation element fixed in the population in most simulations when the chance of infection was 50% or 90%. Networks featuring pleiotropic downregulation fixed in the population when the chance of infection was 50%, while slowly evolving pleiotropic proteins had a high fixation rate at the intermediate infection risk. In evolved competition, networks featuring pleiotropic downregulation at high infection risk and those featuring fixed upregulation at intermediate infection risk fixed in most simulations. All results are shown in Fig C in S1 Text. Collectively, these results suggest that pleiotropic networks are capable of outcompeting non-pleiotropic ones at intermediate and high infection risks over evolutionary time.

Inducible immunity increases fitness relative to constitutive immunity, but is much rarer
The results that we obtained from this model generally favor the evolution of constitutive immunity as predominantly inducible immunity was rare in most scenarios, but it is unclear if this was due to high fitness imparted to hosts by constitutive immunity or the evolutionary feasibility of producing an inducible immune response. The only populations in this study that consistently produced predominantly inducible immune responses were the pleiotropically downregulated hosts at higher infection risk levels (Fig 2A). When comparing the absolute correlation coefficient calculated between the non-pleiotropic and pleiotropic networks. B) Heatmap of the magnitude of maximum immune response attained during infection vs proportion of immune response that is induced by parasites. The y-axis shows the peak of immune effector abundance achieved during infection range [0,1]. The x-axis shows the proportion of the peak response that was generated following infection range [0,1] Darker colors indicate more individuals expressing the magnitude and response combination.
https://doi.org/10.1371/journal.pcbi.1010445.g002 The y-axis shows the mean absolute difference in effector levels between intact networks and single signaling protein knockout networks. Effector levels were measured during infection by a parasite that could not manipulate host immune signaling. Blue dots correspond to pleiotropic signaling protein knock outs, green dots correspond to knocking out of a single non-pleiotropic protein. All plots on the same row have the same chance of infection (10%, 50%, or 90% descending), and plots in the same column have the same implementation of pleiotropy. There are no pleiotropic nodes in non-pleiotropic networks (leftmost column), so nodes were just chosen at random twice. The networks used in this analysis were the most common network at the end of the simulation from which they originated. Asterisks denote significant differences between pleiotropic and non-pleiotropic knock outs.
https://doi.org/10.1371/journal.pcbi.1010445.g003 The immune response of the winning population of competitive scenarios was almost always more inducible than the immune response mounted by the losing population. Immune response probability density function of pleiotropic winners (blue) and the last non-pleiotropic network (green) in the simulation. All plots on the same row have the same chance of infection (10%, 50%, or 90% descending), and plots in the same column have the same implementation of pleiotropy. The x-axis shows the percent of the response that is induced by parasites, with the left-hand side being 0% of response induced, to 100% induced responses on the right. The y-axis corresponds to the relative likelihood of finding an immune response in the specified population that is X% induced. 10 of 12 scenarios show pleiotropic winners being more inducible in their immune responses (as determined by right shifted density peaks) than their non-pleiotropic competitors. Each plot shows the results for competition after 250 generations of evolution. https://doi.org/10.1371/journal.pcbi.1010445.g004

PLOS COMPUTATIONAL BIOLOGY
fitness of our populations we see that pleiotropic downregulation exceeded the fitness of nonpleiotropic hosts at higher infection risk levels, while all other network types were approximately equally fit (Fig J in S1 Text). This suggests that hosts expressing inducible immune responses are more fit than constitutive hosts. To investigate whether this absolute fitness advantage translated to a competitive advantage we looked at competitive simulations between pleiotropic and non-pleiotropic populations focusing on the most prevalent network from the winning population and the last network from the losing population (Figs 4 and D-I in S1 Text). First comparing pleiotropic winners to non-pleiotropic losers, we saw that in 10 of 12 scenarios inducible immunity was more common in the pleiotropic winners than in the nonpleiotropic losers (Fig 4). Non-pleiotropic winners were similarly more inducible than pleiotropic losers (Fig F in S1 Text). No pattern emerged when comparing pleiotropic winners and non-pleiotropic winners or pleiotropic losers and non-pleiotropic losers (Figs H and I in S1 Text). These results suggest that evolved inducible immune responses are more fit regardless of their interaction with pleiotropy, and this is further supported by our findings that the existence of a single host expressing a highly inducible immune response was indicative of many hosts being similarly inducible (Fig N in S1 Text). We propose then that pleiotropy facilitates the evolution of inducible immunity which then bestows higher relative fitness on its host.

Hosts initially express constitutive immunity and transition to inducible immunity over evolutionary time
Examining the transition of response types over time could shed light on the network features that facilitate the evolution of inducible immune responses. When looking at the immune responses and magnitude of response for each population in the first 50 generations of co-evolution we saw that in most scenarios populations rapidly converged on a single equilibrium over time, differing only in the magnitude of effector activated (Figs O-Q in S1 Text). In the pleiotropically downregulated populations, however, we observed that populations branched toward multiple alternative states after converging on an early equilibrium. This difference was especially pronounced when the risk of infection was 90% (Fig Q in S1 Text).
The degree of inducibility and magnitude of immune effector activated by hosts in the final generation did not seem to be restricted by the initial ancestor of a host. When looking at the lineages of hosts present in the last generation of a simulation, we found that hosts in inducible populations often shared a progenitor with hosts in constitutive or mixed-strategy populations (Fig R in S1 Text). While the complexity of our model prevents us from analytically solving for the existence of evolutionarily stable strategies, this behavior provides evidence for the presence of immune response strategies that are, at least temporarily, co-stable.

Discussion
Our model of signaling network evolution allowed us to investigate the changes in evolutionary trajectories and endpoints that are associated with pleiotropic signaling proteins. Our findings suggest that pleiotropy in immune networks can be beneficial to organismal fitness, both by speeding the development of highly fit immune response dynamics and encouraging the exploration of phenotypic space by easing the transition from local fitness peaks to global ones. The effects of pleiotropy on organismal immune responses depend both on the specific action of the pleiotropic protein and on the prevalence of parasitic antagonists in the environment. These factors contribute to an evolutionary landscape where peak immune effector levels depend on parasite abundance, and the mechanisms of achieving those peaks are heavily influenced by pleiotropic signaling proteins.
While the networks generated using this model cannot be explicitly mapped to any specific biological signaling pathway, there are important similarities between simulated host immune networks and those found in nature (Fig S in S1 Text). In both cases, signaling networks that make use of feedback loops between proteins and self-regulatory behavior is common [30,31]. The evolved networks were sparsely connected, a feature often associated with biological signaling networks [32], with the mean connectedness of a population rarely exceeding 50% ( Fig  L in S1 Text). Some explicit parallels can be drawn to Caenorhabditis elegans immune defense, where immune action is carried out by parallel signaling pathways that communicate but are capable of independent action and do not share a single mechanism of control [33], a feature observed in our populations as well (Fig S in S1 Text). We observed that our evolved immune networks mimic biological ones not only in structure but also in function; constitutively immune hosts have lower peak immune effector abundances than hosts with inducible immune systems. This finding is consistent with expected stimulus responses in multilayered regulatory systems [34] and has been observed previously in explicit models of immune signaling networks [35]. Our proposal that pleiotropic proteins can promote the evolution of complex traits, such as inducibility, is rooted in the understanding that extant biological processes are often co-opted when adaptive evolution gives rise to novel traits [22,24,36].
The hosts in our model tend to evolve constitutive immune responses even when the chance of infection is low, a condition that has previously been hypothesized to favor inducible immunity [9,10,37]. We believe this difference arises because we are not solely assessing the relative fitness of these immune strategies, but also their evolvability. The low abundance of parasites simply does not present enough of an evolutionary pressure for hosts to develop inducible immune networks. As parasite abundance increases, constitutive hosts carefully balance their effector abundance against the rate of parasite replication to prevent growth, but rarely entirely remove parasites, in line with behavior predicted by Shudo and Iwasa [38]. We see evidence of this balancing act in feedback loops between proteins that upregulate the effector and the effector itself, as well as effectors upregulating proteins that suppress their own activity (Fig S in S1 Text). Constitutive immunity then represents a local fitness maximum that is easy to attain, especially for networks early in the evolutionary process (Figs O-Q in S1 Text).
Inducible immunity is generally a rare response in the context of our results, which is in line with previous work that modeled signaling networks explicitly [11] as well as those that developed analytical models of optimal trait usage [5,6,37]. We might naively expect, then, that constitutive immunity would be more fit than inducible immunity. However, our results suggest that the opposite is true-populations that produced predominantly inducible responses were significantly more fit than their constitutive counterparts at higher infection risk levels (Fig J in S1 Text). When looking at population immune responses in the first 50 generations of a simulation, we see that hosts in the first generation predominantly mount constitutive responses across all conditions studied and that the first predominantly induced hosts do not appear until later generations (Figs P and Q in S1 Text). The scarcity of inducible immunity and the extended evolutionary time necessary for highly inducible immune responses to arise reinforces the notion that inducible immune responses are evolutionarily complex to deploy. We believe that fixed downregulatory pleiotropy results in the evolution of inducible responses more often than for fixed upregulatory pleiotropy because while the latter results in overinvestment in immunity it still provides protection against infection. Uncontrolled downregulation, on the other hand, silences the immune response and leads to infection-induced host death. The resulting selective pressure favors the evolution of upregulatory architecture to compensate, paving the way for inducible immunity.
It has previously been speculated that part of the evolutionary challenge of evolving inducible immune responses may be that they are more susceptible to disruption via parasites than constitutive immune responses [11]. This thinking is supported by our data, especially when parasite manipulation targets the pleiotropic signaling protein (Fig 3). Based on these findings we suspect that the early evolutionary steps leading to inducible networks can be susceptible to manipulation via parasites, potentially negating early fitness gains over constitutive hosts. If this were the case, then once the hurdle of manipulation is overcome (e.g., with signaling redundancy) that fragility is reduced, and inducible hosts are left with a highly fit immune response that requires energetic investment only in times of infection. Another possible explanation is the accumulation of nodes over time that decrease parasite fitness if manipulated, forcing the parasites to avoid architecture that is critical to the induced immune response. As host networks grow, moreover, the proportion of the network affected by any given parasite manipulation shrinks, so networks that grow to a critical size could minimize the probability of critical parasitic manipulation during an induced immune response.
Immune effectors across the tree of life are incredibly diverse in their form, function, and targets, making it impossible to generalize their behavior in a tractable model. Therefore, we made the following simplifying assumptions with regards to the immune effectors implemented in our model: effectors molecules produced by hosts were perfectly effective at removing parasites (i.e., an effector abundance of 50% would remove 50% of the current parasite population in a single time step), and that effectors acted immediately to kill parasites. The first of these assumptions could be violated by the non-linear relationship between density and parasite killing activity exhibited by some effectors, including AMPs [39], while the second ignores the process of activation many immune effectors undergo prior to fighting infections, a critical step necessary for preventing autoimmunity [40]. One potential continuation of this work would be to incorporate the relationships between inducibility, non-linearity, and pleiotropy by modifying parasite killing based on effector concentration. To better capture the role of activation, future studies could incorporate a timestep delay so that the parasite population at time t are killed by effectors at time t-n, where n is the desired delay. These types of modifications would facilitate a more accurate evolutionary picture of particular signaling pathways, although at the cost of generalizability across pathways, effectors, and species.
In conclusion, we have identified distinct changes in the trajectory of signaling network evolution associated with the inclusion of pleiotropic signaling proteins. Fixed random pleiotropy and evolutionary rate constraints on the pleiotropic protein did not result in significantly different evolved networks when compared to the non-pleiotropic control. Fixed upregulatory and fixed downregulatory pleiotropy altered initial and terminal network dynamics (Figs 2 and O-Q in S1 Text), network size, connectivity, and the number of distinct paths from the detector to the effector (Figs K-M in S1 Text), all while maintaining mean population fitness that was approximately equal to or greater than non-pleiotropic hosts (Fig J in S1 Text). With these findings we have provided some of the first evidence for the wide-ranging evolutionary effects of pleiotropic proteins on the signaling networks they are a part of, highlighting the importance of further empirical investigation into the benefits, tradeoffs, and evolutionary consequences of pleiotropy in the immune system and across the genome. This work advocates for using a broad perspective when studying known pleiotropic proteins and genes, as their full evolutionary effects may only be observed at the scale of signaling networks.

Materials and methods
The model of signaling network evolution that we based this work on was originally designed to model the evolution of generic signaling networks [41], and was subsequently modified to study immune signaling network robustness [42] and the evolvability of immune responses during coevolution [11]. Here we have revised the model to include pleiotropic signaling proteins. We designed the following implementations of pleiotropy: 1) random connections between the pleiotropic protein and other proteins in the network that are fixed across evolutionary time; 2&3) a single connection from the pleiotropic protein that up-or downregulates the effector as a side effect of its function in development; 4) pleiotropic proteins that may evolve but undergo reduced mutation rates relative to other network components, capturing reduced evolutionary rates associated with purifying selection . Fig 1 provides a diagrammatic representation of these restrictions compared to the non-pleiotropic case.
We used two broad classes of simulation to study the effects of pleiotropy on immune evolution: co-evolution and competition. In co-evolution simulations, a population of hosts evolved for 500 generations with a population of parasites. The host immune networks in these simulations were either non-pleiotropic or all hosts in the population had the same pleiotropic constraint as defined above. Each simulation had 500 hosts, with each host initially defined by a randomly generated immune network. Parasite population size was determined by the chance of infection in each simulation, and each parasite possessed a single connection to a signaling protein that could disrupt host immune signaling. We used these simulations to study immune networks and their dynamics when evolving under pleiotropic constraint.
Competitive simulations were broken into two phases: independent evolution and competition. During independent evolution a non-pleiotropic population (n = 250) and a pleiotropic population (n = 250) were allowed to co-evolve with parasite populations without interacting with the other population of hosts. After 250 generations the simulation entered competition, combining the host populations along with their parasites. The competition ended when one population died out entirely or 1000 generations had passed with no winner (draw). These simulations allowed us to evaluate the viability of pleiotropy in a population that is not uniformly facing the same potential fitness deficits. In all cases, for each implementation of pleiotropy and chance of infection, we conducted 100 simulations.
Each host network initially contains a single detector, three signaling proteins, and a single effector. Over the course of a simulation, mutations during reproduction duplicated or deleted signaling proteins and deleted, added, or altered regulatory interactions between proteins in the network. Hosts remained restricted to a single detector and a single effector, and at no point were detectors and effectors allowed a direct connection. Host fitness was determined as a function of immune effector abundance pre-and post-infection, cumulative parasite load, and network size (see Eq 3). Parasite fitness was strictly based on cumulative parasite load during infection, a proxy for transmission potential.

Simulation Framework
Evolutionary simulations were carried out in the Julia programming language (v 1.6) and can be broken into the following discrete events.
1. Initialization: a population of hosts is generated at random. All hosts start with a detector, three signaling proteins, and an effector. Each possible connection in the network is given a 50% chance of occurring, and for each connection that does occur a random number is selected on the interval [-1,1] to decide the regulatory behavior of that connection. The only constraint on initial network structure is that the detector protein cannot directly connect to the effector protein.
A population of parasites is also generated to infect 10%, 50%, or 90% of hosts, depending on the chance of infection for the simulation. Parasites are designated by a target (restricted to signaling proteins) and regulatory behavior on that target (restricted to the interval [-1,1]).

Healthy Equilibrium:
Each protein P in an immune network has total concentration equal to 1, with the active portion denoted P i � and inactive portion denoted P i , ([P i � ] + [P i ] = 1). When determining the effects of protein P on other proteins in the network only the active portion is considered. During infection, changes in parasite abundance are calculated as though it was another protein in the network. The change in P � i is determined by the regulatory action on P i defined: Where k i,j are the upregulatory coefficients from protein P j to protein P i , I i,j are the downregulatory coefficients from protein P j to protein P i . The term Protein Use describes the inactivation of some P � i owing to the action of P i on other proteins in the network. Starting from initial concentrations of P � = .5 for all proteins in the network, P � for each protein at time step t+1 is calculated numerically using Eq [1] until the effector protein reaches equilibrium.

Fitness Calculation:
Using data from the Healthy Equilibrium and Infected Equilibrium phases, fitness is calculated using the following equation: With ( ½P � eff � pre refers to the amount of active immune effector prior to infection, capturing the cost of constitutive investment in immunity. ½P � eff � post refers to the amount of active immune effector at the end of the simulation and would penalize an inducible response if it was not adequately regulated. The variable V is a damage coefficient that is akin to parasite virulence, Area is the area under the parasite infection curve normalized to 1, and ProtCost controls for a loss of fitness associated with the number of proteins in the network. The network size cost used the size of the Toll pathway as a neutral benchmark [43] and the increase in cost with increasing size reflects the fitness deficits associated with maintaining and producing large amounts of proteins [25]. In this model there are two main costs of immunity to fitness: energetic and immunopathological. Over-investment in immunity, either constitutively ð½P � eff � pre Þ or through poorly regulated inducible responses (as captured by continued investment after parasite clearance, ½P � eff � post ), will therefore trade off with other life history traits and host survival [1,44]. The virulence term contains fitness costs due to parasite-induced pathology. By integrating the collective costs of energetic investment, immunopathology, and virulence into a single fitness function, we can evaluate the fitness effects of immune dynamics without forcing an unnatural binning of responses into strictly constitutive or inducible immunity.
Parasite fitness was derived from the normalized area of the parasite infection, which is representative of the total number of parasites produced over the course of the infection. 6. Host and Parasite Death: in each generation up to 30% of each population dies, with hosts that succumbed to infection taking priority in the death process. In the case where greater than 30% of either population died due to poor fitness or uncontrolled fitness, the excess deaths were added back to the healthy population following random selection. All individuals that survived a generation replaced themselves (as well as any other offspring they produced in Step 7) in the next generation. Population size and deaths were capped as a concession to the low fitness of initial randomly generated networks and computational expenses.
If less than 30% of the population died due to parasite burden, surviving individuals were selected to die in a fitness weighted manner. At random an individual was selected from the population and its chance of dying was inversely proportional to its relative fitness against the population. Ex: a host that was in the 75 th percentile of fitness has a 25% chance of dying. This process continued until either all surviving hosts had been evaluated once, or until 30% of the population had died. 7. Host Reproduction: survivors are picked in a fitness weighted manner to reproduce, with a host's chance to reproduce being inversely proportional to its relative fitness against the population. A single host could produce multiple offspring in a reproductive stage. The host population was completely replenished in each reproductive stage (keeping population size constant across generations). Reproduction results in the creation of a direct copy of the parent or, rarely, a mutated copy (host mutation rate: 5e-3). For hosts, a mutation could be one of the following modifications to the parent network: 1) Add a protein-protein interaction between two randomly selected proteins with a random regulatory behavior (relative probability = .25) 2) Delete a protein-protein interaction (relative probability = .25) 3) Alter regulatory coefficient (relative probability = .3) 4) Delete a protein (relative probability = .1) 5) Duplicate a protein (relative probability = .1) � � Duplication is the only mutation that can act on the pleiotropic signaling protein. The duplicated pleiotropic protein is treated as a non-pleiotropic signaling protein.

Parasite Reproduction:
surviving parasites are picked to reproduce in a purely fitnessbased manner, with highly fit parasites producing more offspring than their lower fitness peers. This corresponds to higher cumulative parasite load leading to more offspring in the following generation. Parasites with a cumulative load between 0 and .33 produced one offspring, those with a cumulative load between .34 and .66 produced two offspring, and finally a cumulative load between .67 and 1 lead to 3 offspring. Parasites reproduced until the population is completely replenished. Parasites reproduced by way of direct copy of the parent, or rarely a mutated copy was produced (mutation rate of 1e-2 for parasites). Parasite mutations are defined as follows: 1) randomly change the target signaling protein (relative probability = .5) 2) change the regulatory behavior on the targeted signaling protein to a new value in [-1,1] (relative probability = .5) After step 8, the simulation repeats starting at step 2 using the hosts and parasites populations produced in step 7 & 8. Each of these cycles from step 2-8 is termed a generation and each simulation runs for 500 generations. For each combination of infection rate and pleiotropic constraint we ran 100 simulations. This number of simulations was chosen to balance computational time against reproducibility.

Competitive simulations
We devised competitive simulations to determine the relative fitness differences between pleiotropic and non-pleiotropic hosts in this model. We conducted 100 competitive simulations for each implementation of pleiotropy and chance of infection pairing, and these simulations had the following changes from the non-competitive case described above: 1. 250 pleiotropic hosts and 250 non-pleiotropic hosts either immediately entered competition (unevolved competition) or were allowed 250 generations to evolve independently (evolved competition), at which point their populations were merged and they entered competition.
2. Competitive simulations proceeded until one of the host populations was extinct, resulting in a victory for the extant population, or until 1000 generations had passed with no winner, resulting in a draw.

Model Assumptions
One of the critical assumptions we made in the construction of this model was that pleiotropic genes are immutable except in the case of slowed evolution. Pleiotropic genes in nature are not necessarily immutable, but they tend to evolve at a slower rate than non-pleiotropic genes [15,18,19]. We decided that making the pleiotropic proteins immutable was the best way to capture the disparity in the rate of evolution between pleiotropic and non-pleiotropic proteins in the relatively short generation time we used as a concession to computational limits, although we relax this assumption in the case of the slow evolution condition. We also decided to preserve pleiotropic immutability throughout a simulation, even following a duplication event, despite evidence suggesting that pleiotropic proteins can lose their pleiotropic nature (thus their immutability) through subfunctionalization following duplications [45]. Due to the suspected scarcity of subfunctionalization events [46,47] and the complexity that appropriately implementing subfunctionalization would add to the model, we felt justified in a compromised duplication behavior that preserved the immutable nature of the original protein but allowed the paralog to accumulate mutations.

Data Analysis
Immune response probability density. To study the immune responses mounted by host populations, we first calculated the percentage of each host's maximal immune response that was activated by infection. For each host this generates a value describing the degree of inducibility of their immune response, with 0% indicating no induced response to parasites, and 100% indicating an entirely induced response. We then approximated the probability density function for this data using kernel density estimation. This immune response probability density conveys the likelihood that a host in a population would have a specific percentage of their immune response induced by parasites and these values were normalized to one to ease comparisons between populations. All infected hosts in the final generation of each simulation were used in the calculation of immune response probability density functions.
Correlation. We calculated the Pearson correlation coefficient to aid in rigorous comparisons between pleiotropic and non-pleiotropic host immune response densities. For each combination of pleiotropic implementation and chance of infection, we calculated the Pearson correlation coefficient between pleiotropic immune response density and non-pleiotropic immune response density at the same chance of infection. Incidents of poor correlation (corr�|.2|) suggest a significant difference between the immune responses mounted by the pleiotropic and non-pleiotropic populations.
Immune effector abundance vs immune response type. To visualize the relationship between immune response type and peak immune effector abundance we calculated a twodimensional probability density function, where the x axis was the proportion of the response induced by parasite and the y axis was the maximum amount of immune effector activated. These probability density functions were calculated for each implementation of pleiotropy and chance of infection pairing using kernel density estimation. All hosts infected in the final generation of each simulation were used to generate these plots. The two-dimensional probability density functions for each pairing of pleiotropic implementation and chance of infection were then added together to produce Fig 2B. Determining the effect of signaling protein knockouts. To determine if pleiotropic hosts were more susceptible to manipulation than non-pleiotropic hosts, we calculated the mean absolute difference in active effector levels between intact hosts and hosts with a single signaling protein knockout (the protein was removed from the network). Hosts were infected with a non-disrupting parasite (a parasite that could not interfere with host signaling proteins) for twenty time-steps and the effector levels for the intact network and the knockout for each signaling protein were measured. The mean of the absolute difference in effector levels at each time step between the intact and knockout networks was calculated and is used as a metric of the networks reliance on a specific signaling protein to produce the evolved response. Significant differences between the mean absolute difference in effector level following knockout of the pleiotropic signaling protein was compared that of non-pleiotropic knockout using twotailed homoscedastic t-tests with Bonferroni correction.
Network features. Network connectivity was calculated by dividing the number of protein-protein interactions in a network and dividing that number by the number of possible connections that network could possess. To determine if pleiotropy altered the number of proteins necessary to mount an immune response, we measured network size by counting the number of proteins present in the most common immune network from the end of each simulation for a given implementation of pleiotropy and chance of infection. Finally, as a measure of robustness, we calculated the number of distinct paths from the detector to the effector. A distinct path was one that did not share a signaling protein with any other path connecting the detector and effector, and previous empirical work has proposed this to be a mechanism promoting robustness to parasitic interference in immune network signaling [33].
Supporting information S1 Text. Table A. Names, values, and description for variables and parameters used in the simulation. Fig B. Diagram of infection period, end states of infection, and example of how key findings were drawn from infection data. A) the host is infected, and the infection dynamics are calculated as described in the Methods Simulation Framework step 4. Infection ends in one of 3 ways: B) the parasite is managed, but not killed before the 20 timesteps have passed, C) the parasite is killed before the 20 step limit is reached, D) the parasite goes unmanaged and kills the host after 20 time steps have passed. E) shows an example of how the data used to generate immune response density plots were collected. Following the conclusion of the infection, the difference between initial effector abundance and maximum effector abundance was determined for each infected host. The percent of the maximal abundance that was induced by parasites was then calculated for each host and used to generate an immune response probability density plot for the population. Created with BioRender.com. Fig C. Pleiotropic hosts can outcompete non-pleiotropic hosts. Results of competition simulations. Rows correspond to infection percentages and columns correspond to the pleiotropy type for a set of competitions. Unevolved competitions are those that had non-pleiotropic and pleiotropic organism enter competition immediately. Evolved are those that took place after 250 generations of adaptation in isolated populations. Fig D-I. Winners of competition simulations are consistently more inducible than corresponding losers, but pleiotropic and non-pleiotropic hosts are similarly inducible when matched for winning and losing. The immune response density plots of winners and losers of competitive simulations after 250 generations of adaptation. Comparisons presented are a) pleiotropic winners vs. non-pleiotropic losers, b) pleiotropic winners vs. pleiotropic losers, c) non-pleiotropic winners vs. pleiotropic losers, d) non-pleiotropic winners vs. non-pleiotropic losers, e) pleiotropic winners vs. non-pleiotropic winners, f) pleiotropic losers vs. non-pleiotropic losers. Fig J. Down regulatory pleiotropy results in hosts for whom fitness that equals or exceeds non-pleiotropic hosts. Plots of average host fitness through 500 generations when the chance of infection was 10%(left), 50% (middle), or 90% (right). Each panel shows host or parasite fitness from unconstrained (solid line), Fixed Random (squares), Fixed Up (triangles), Fixed Down (diamonds), and 100x slower evolution (circles) simulations. Average host fitness was calculated using hosts that were and were not infected for each generation. Note that the y axis changes scale in the second and third panel to because overall host fitness decreased as the chance of infection increased. Fig K. Pleiotropy can significantly alter the size of host signaling networks. Size (number of proteins) of the most common network from each run of a scenario with median lines presented in black. All networks start with 5 proteins, and at the 10% infection rate, no pleiotropic conditions diverge significantly from this point. At the 50% infection level, only the fixed downregulation pleiotropic condition has an increased median network size. Conversely, at the 90% infection level, only fixed upregulation did not increase in median network size. Asterisks indicate a significant difference from the non-pleiotropic scenario in each row. Fig L. Pleiotropy can significantly reduce signaling network connectivity. Percentage of total potential connections deployed by the most common network at the end of each simulation. Each row is an infection risk (10, 50, or 90%) and each column is a type of pleiotropy. On average, half of all connections are used initially. Asterisks indicate a significant difference from the non-pleiotropic scenario in each row. Fig M. Pleiotropic hosts can develop significantly more distinct paths through a network than non-pleiotropic hosts. Number of distinct paths from the detector to the effector in the most common network following a simulation. The distinct paths through a network are the set of paths that share no signaling proteins with the other paths in the set. The fixed downregulation conditions deploy a higher number of distinct paths. Importantly, distinct paths connect the detector to the effector in a manner that is partially insulated from other paths through the network, increasing robustness. Asterisks indicate a significant difference from the non-pleiotropic scenario in each row. Fig N. A single highly inducible host in a population is often indicative of many hosts that are at least that inducible. Plots show the average percentage of a population that is at least as inducible as the level indicated on the x-axis. X-axis is the inducibility threshold and the Y-axis is the average proportion of the population that meets or exceeds that threshold given that at least one host meets or exceeds it. Columns correspond to the implementation of pleiotropy, rows correspond to the percent of the population that is infected in each generation (10, 50, 90 percent) Example from the Fixed Down column, 90% infection row: Looking at the bar labeled .9, this indicates that in all simulations with at least one host where > = 90% of the maximal immune response was induced by parasites a significant majority of hosts (~90%) where also expressing immune responses that were > = 90% inducible. Fig O. When the chance of infection is low, the evolutionary trajectory of hosts does not depend on their pleiotropic status. Magnitude of immune response by the proportion of response that is induced in the initial 50 generations of an evolutionary simulation where the chance of infection was 10%. Darker colors indicate more common combinations of magnitude of immune responses and proportion of response induced by parasites. Fig P. When the chance of infection is moderate, downregulatory pleiotropy leads hosts to novel evolutionary trajectories. Magnitude of immune response by the proportion of response that is induced in the initial 50 generations of an evolutionary simulation where the chance of infection was 50%. Darker colors indicate more common combinations of magnitude of immune responses and proportion of response induced by parasites. Fig Q. When the chance of infection is high, downregulatory and upregulatory pleiotropy led hosts to novel evolutionary trajectories. Magnitude of immune response by the proportion of response that is induced in the initial 50 generations of an evolutionary simulation where the chance of infection was 90%. Darker colors indicate more common combinations of magnitude of immune responses and proportion of response induced by parasites. Fig R. As the chance of infection increases, the proportion of host lineages that contain constitutive and inducible hosts increases. The proportion of runs where hosts that descended from the same initial host ended up with immune systems with different response dynamics (i.e., a split or bifurcated lineage) in the final generation of the simulation. We refer to hosts that share an ancestor but do not share immune response dynamics as being a part of a split lineage. Fig S. The pleiotropic nature of a host shapes end state signaling networks. The average host network generated in each pleiotropic constraint and infection level pairing. These average networks were generated using the most common networks from the end of each simulation at a given pairing. Saturation for the connections between proteins is scaled based on the most common connection across all networks at the given constraint and infection level. Arrows denote the direction of the connection, blue connections are down regulatory, red are upregulatory.